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ABSTRACT 

We present a numerical method to estimate the lensing parameters and the Hubble 
constant Hq from quadruply imaged gravitational lens systems. The lens galaxy is 
modeled using both separable deflection potentials and constant mass - to - light ratio 
profiles, while possible external perturbations have been taken into account introduc- 
ing an external shear. The model parameters are recovered inverting the lens and the 
time delay ratio equations and imposing a set of physically motivated selection criteria. 
We investigate correlations among the model parameters and the Hubble constant. Fi- 
nally, we apply the codes to the real lensed quasars PG 1115H-080 and RX J0911-I-0551, 
and combine the results from these two systems to get Hq — 56±23 km s^^ Mpc^^. 
In addition, we are able to fit to the single systems a general elliptical potential with 
a non fixed angular part, and then we model the two lens systems with the same 
potential and a shared Hq: in this last case we estimate Hq = 491^]^ Km s~^ Mpc~^. 
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1 INTRODUCTION 

Gravitational lensing is one of the main tools to obtain in- 
formation about the structure and evolution of the universe. 
In particular, time delay measurements are a recent primary 
distances indicator, furnishing a new method to estimate 
the Hubble constant Ho, which determines the present ex- 
pansion rate of the universe (see, for instance, (Narayan & 
Bartelmann 1998; Kochanek & Schechter 2003)). 

Actually, in 1964, Refsdal proposed to estimate Hq 
(Refsdal 1964a; Refsdal 1964b; Refsdal 1966) from multi- 
ply imaged QSOs by the measurements of the delays in the 
arrive time between light rays coming from the different im- 
ages, which follow different optical paths. It is not difficult to 
show that the time delay between two images due to a grav- 
itational lens can be factorized in two pieces: the first one 
depends on cosmological parameters and is inversely propor- 
tional to Hq, while the second one is determined by the lens 
model only. Thus, having measured time delays among im- 
ages of a lensed QSO, once we fix the cosmological parame- 
ters, we can obtain a direct estimate of -ffo provided that the 
lens model has been recovered from the lensing constraints, 
or is known in an other independent way. Nowadays, there 
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are more than sixty multiple image systems (CASTLES web 
page), but only about ten of them have measured time de- 
lays. However, this number is increasing day by day, and 
in the future it will be possible to measure time delays for 
many other systems. 

It is well known hat the most significant uncertainty af- 
fecting the estimate of Hq with the Refsdal method is only 
related to the mass model used. In the usual approach the 
model parameters are recovered by fitting some parametric 
models to the available constraints through minimization 
techniques. Instead, other authors (see (Williams & Saha 
2000)) carried out the "pixellated lens" modelling, that de- 
scribes the mass distribution by a large number of discrete 
pixels with arbitrary densities, so determining the Hubble 
constant by means of a set of physical motivated constraints. 
A compromise between these two approaches consists in the 
numerical solution of a set of non linear equations: in a pre- 
vious paper ((Cardone et al. 2002), hereafter CCRP02), for 
instance, some of us applied a semianalytical technique to 
fit the separable potential of the form if; = r"F{6) and were 
able to develop an algorithm to estimate Hq without the 
need to give an explicit expression for the shape function 

Here, we further improve the general method in (Car- 
done et al. 2001) (hereafter CCRPOl) and CCRP02, imple- 
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meriting a set of Mathematica codes to shape quadruple lens 
systems and obtain a better estimate of the Hubble constant. 
Such a method in fact allows us to manage a still wider class 
of lens models and to obtain useful information about lens 
parameters and Hq, eliminating possible not physical solu- 
tions. 

Actually, we consider both separable models, specifying 
the angular part F{9), and models with constant mass-to- 
light ratio. As usually in literature, we develop the poten- 
tial of an external object contributing to the Icrising phe- 
nomenon to second order, and hence its effect on the total 
deflection potential translates into an external shear term 
(Schneider, Ehlcrs & Falco 1992). Wc constrain the lens 
models using the image positions with respect to the lens 
centre and the time delay ratios, which allow us to highly 
constrain the used models, even if the equations to be solved 
become more complex. Once the results from all the possi- 
ble models are collected together, such a procedure offers 
the advantage of giving an estimate of the Hubble constant 
which is in a sense marginalized with respect to the lens 
models. 

Moreover, our numerical method allows to pursue a new 
kind of approach to the lens modelling, that considers at the 
same time several lens systems. Actually, following (Saha & 
Williams 2004), we try to model two lensed systems with the 
elliptical potential tp = r"F{9) and an external shear, based 
on a strong hypothesis we make ab initio: we suppose that 
the lens systems has a shared Hq. In this way we are able 
to create a two-system model that can give a more complete 
estimate of Hq. 

The outline of the paper is as follows. In Sect. 2 we write 
the lens equations and the time delay function in polar coor- 
dinates. A careful presentation of the models wc use is given 
in Sect. 3, while Sect. 4 is devoted to the presentation of the 
numerical method used to 'fit' single lens models to systems. 
We also discuss the constraints used to select physically mo- 
tivated solutions and the statistical approach to obtain the 
final estimate of the parameters. Then in Sect. 5 we present 
a new approach to shape two lens systems with a shared 
Ho. In order to verify whether our codes are able or not to 
recover the correct values of parameters, we build simulated 
systems to which we apply the codes as described in Sect. 6, 
where, in addition, we analyze the biases and the uncertain- 
ties in the use of a model respect to another one. In Sect. 7, 
we use the simulated systems to investigate the existence of 
degeneracies among the model parameters, recovering some 
interesting scaling laws. Sect. 8 is dedicated to the appli- 
cation of our procedure to two real quadruple systems, PG 
1115-1-080 and RX J0911-I-0551: we obtain an estimate of 
the lensing parameters and the Hubble constant, also taking 
into account the contribution of changing the cosmological 
parameters to the uncertainty on Ho. Finally, we present a 
discussion of our results and future improvements in Sect. 
9. 



the rectangular ones through the following coordinate trans- 
formation: 

a; = — rsin^, y = r cos 9. (1) 

Time delay function, i.e. time delay of a generic path from 
source to observer, is given by (Blandford & Narayan 1986; 
Zhao & Pronk 2001): 



At = h noo 



L J.2 _ j,^^ cos{9 - 



T 2 s 



i,{r,e) 



(2) 



where {r,9) determines the impact position of the generic 
path on the lens plane (with r measured in arcsec), {rs,9s) 
is the unknown source position, and tp{r,9) is the adimcn- 
sional lensing potential. Then, h is the Hubble constant Ho 
in units of 100 km s~^ Mpc~^, while noo is a typical time 
delay linked to cosmological parameters, and defined as: 



TlOO = 



DgPs 1 + Zd 

Dds c 



(3) 



where Dd, Ds and Dds arc observer- lens, observer- source, 
and lens - source angular diameter distances, calculated with 
Ho = 100 km s~^ Mpc"'^, and Zd is the lens redshift. This 
factor contains all the cosmological information, since the 
distance depend on the other cosmological parameters: i.e. 
the density matter parameter Sim, the "dark energy" con- 
tribution Qx (see (Peebles & Ratra 2002; Caldwell, Dave & 
Steinhardt 1998; Demianski et al. 2003)) and the smoothness 
parameters a introduced in (Dyer & Roeder 1972; Dyer & 
Roeder 1973; Dyer & Roeder 1974) to take into account the 
inliomogcneities of the universe. If i and j are two images, 
the time delay between them is AUj = AU — Atj . 

According to Fermat principle, the images lie at the 
critical points of At, so one can obtain lens equations mini- 
mizing it. We get: 



r — re cos{9 — 9s) = 



dr ' 



(4) 



rssin{9-9s) = \^. (5) 

We shall use a lensing potential formed by the sum of two 
terms: 



V'(r, 9) = ipier^s (r, 9) + ipext{r, 9), 



(6) 



where ipuns {r, 9) is the contribution of the lens galaxy, and 
^Ext{r, 6) is the external perturbation. The first term is con- 
nected with the mass distribution of the galaxy through the 

Poisson equation: 



y^i^le-na{r,e) = 2K(r, 



(7) 



being K(r,9) the convergence, i.e., the adimensional surface 
mass density, defined as: 



K{r,e) 



(8) 



2 LENS EQUATION AND TIME DELAY 

Let us choose a rectangular system {x,y) centered on the 
lens galaxy and with axes pointing towards West and North, 
respectively. Let (r, 6) be the polar coordinates, being 9 the 
position angle measured from North to East, connected to 



where T,crit = 4„oj^o — • The second term describes the 
effects on the lensing system of the environment, i.e., of the 
cluster of galaxies which the lens galaxy belongs to, or an 
external group of galaxies. We describe this contribution 
developing the lensing potential of the external perturbation 
to the second order: 
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Table 1. Separable deflection potentials. 



Model 




Model 1 


cx 






Model 2 


b^ysin{e - dqY + 7-2 cos(6» - e<,)2 r 






Model 3 


^sin{e - 65)2 + g-2 cos(6» - OqY r« 



ipext{r,9) = ^l)shear{r,9) = - -^r^ cos2{e - 6^) , (9) 

where 7 is the external shear and 6~f is the shear angle, ori- 
ented from North to East and pointing towards the external 
perturbation. 



3 LENS MODELS 

The majority of the lens galaxies are early-type galaxies, 
since lensing selects galaxies by mass and the fact that early 
type galaxies are more massive than spirals overwhelms the 
fact that spirals are slightly more numerous. Early type 
galaxies show a wide variety of optical shapes including 
oblate, prolate and triaxial spheroids. Moreover, also their 
mass distribution is not yet fully understood: stellar dy- 
namics and X ray observations all suggest that the early 
type galaxies are dominated by dark matter halos; on the 
other hand, some dynamical studies recently performed us- 
ing planetary nebulae as tracers in the halos of early-typo 
galaxies show evidence of a universal declining velocity dis- 
persion profile, and dynamical models indicate the presence 
of little dark matter within 5iie - implying halos either not 
as massive or not as centrally concentrated as CDM pre- 
dicts (Romanowsky ct al. 2003) . So not only a fundamental 
hypothesis of the CDM paradigm have been remained up 
to now largely unverified - that there should be similarly 
extended, massive daxk halos around ellipticals -, but pre- 
dictions about the detailed halo properties have not been 
testable. Gravitational lensing is actually an unique tool to 
studying elliptical halos, and face the question about the 
nature and the distribution of dark matter in early type 
galaxies, with implications on the estimation of the Hubble 
constant. In order to consider both a luminous-dominated 
and dark-dominated coiiipouciit in the main lens galaxy, we 
consider two different classes of models : separable potentials 
and constant mass - to - light (M/L) mass profiles. 

3.1 Separable potentials 

For the first class, we assign the lensing potential in a simple 
separable form: 

tir,e)=r" F{e,q,eq) , (10) 

where we have explicitly indicated the dependence of the an- 
gular part F{9, q, 9q) on the axial ratio g (0 < g < 1) and the 
position angle 9q that specifies the orientation of the lensing 
galaxy. The potentials (10) axe a generalization of the pseu- 
doisothermal elliptical potentials (PIEP) studied by (Kovner 



1987). This kind of model is in many aspects similar to the 
pseudoisothermal elliptic mass distribution (PIEMD) mod- 
els, such as the singular isothermal ellipsoid (SIE) (Bland- 
ford & Kochanek 1987). However, while PIEMD models can 
represent mass models with any elongation, PIEP models 
represent physically plausible lenses only for some values of 
the axial ratio g, such as g > go, being go a suitable value 
of the axial ratio of the isopotential contours (Kovner 1987; 
Kassiola & Kovner 1995). For q < qo, the elliptical isoden- 
sity profile may in fact be boxy or disky. In particular, when 
dealing with these elliptical models, we adopt the following 
expression for the angular part: 

F oc [sin^(6l - Oq) + q-'^ cos^(6» - Oq)] , (11) 

where 6q is oriented from West to North. The angular part 
F also depends on a strength parameter b: for a singular 
isothermal sphere (a = 1 and F =const) this parameter 
depends on the cosmological parameters and the redshifts 
by means of a distance ratio, and on the central velocity 
dispersion. In our case we can assume a similar dependence 
only for Model 3 (see Table 1). 

In Table 1 we give the expression of for the three 
different models we consider. Model 1 is spherically sym- 
metric, and hence the shape function is F = . On the 
other hand, both Model 2 and 3 are ellipticals, but for Model 
^ we fix a = 1 as for isothermal mass distributions, while a 
is unfixed for Model 3. 

It is worth to note that, since the lensing potential tp 
and the adimensional surface density k are related by a dou- 
ble integration, the ellipticity of the isopotential contours is 
different by that of isodensity ones. Therefore, we have to 
obtain a relation between the axial ratio q of the potential 
and the isodensity axial ratio, named g«. For Model 2, the 
convergence is: 

K{r,d) = ^{coS^ {9- eq)+q^ sin {9 ~eq))~i . (12) 

It is easy to see that k is always positive, i.e., Model 1 is 
always correctly defined, and the analytical relation g« = 
g^ holds. Things get more complicated for Model 3. The 
convergence now turns out to be: 

K{r,9) = ^-^^J^^^[2{-l + q'' +0^)005^ {9- 9q) 
4qa 

-hg^(2(l + q^{-l + a^)) sm''{9 - 9q) + a sin* 2{9 - 9q))] 
X [(cos^(6l - 9q) + q sin^(6> - 6>g))]-^ . (13) 

It is possible to see that, if one of the two factors g^(Q:^— 1)-|-1 
and g^ -f — 1 is negative, then is k < 0. In particular, 
if Q > 1, we always have k > 0, since the two quantities 
above are positive; instead, when a < 1, the convergence 
K is negative if g < \/\ — . For instance, if a = 0.5 and 
g < 0.866, the convergence is negative. The axial ratio g^ 
satisfies the relation 

/g(-l + g2 + Q2)\ 7^ 

9- = 2.1. 22 ' 14 

y — g2 -f 1 -f g2a2 y 

that corresponds to more rapid trends for lower values of a. 
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Table 2. Mass distribution for constant Af/L models. For the de 
Vaucouleurs model see Keeton & Kochanek 1997. 



Model 


K 


Model 4 




Model 5 


1 

2 s^+r^ 



3.2 Constant M/L models 

The second class we consider contains models that describe 
luminosity profiles of elliptical galaxies. We use the dc Vau- 
couleurs (de Vaucouleurs 1948) and Hubble (REF) models 
(respectively Model 4 and 5), whose convergence n is given 
in Table 2. We assume that the mass density is spherically 
symmetric, so that we can avoid any difficulties in solving 
the lens equations for these models. Given the surface den- 
sity K, the deflection angle is easily obtained as (Schneider, 
Ehlers & Falco 1992; Keeton 2001): 



a{r) 



r' dr' K{r'), 



(15) 



and then the deflection potential is evaluated solving the 
equation a = Vip. 



4 SINGLE SYSTEM MODELLING 

Developing the general method used in CCRPOl and 
CCRP02 to model quadruply lensed systems, we numeri- 
cally solve systems of nonlinear equations, imposing on the 
sample of solutions some criteria to select only the physi- 
cal ones, that have to be collected together by means of an 
appropriate statistical analysis. 



4.1 Method 

Here, we consider separable models with an assigned an- 
gular part and more complex mass models for which it is 
necessary to assign the surface density k. In order to take 
into account all of the information that come from a lensing 
event, we make use of time delay ratios, that allow to con- 
strain the lensing parameters quite efficiently. We do not use 
flux ratios since these may be contaminated by microlens- 
ing (Chang & Refsdal 1979; Koopmans & de Bruyn 2000), 
and other effects, such as those due to substructures in the 
lens galaxy (Mao & Schneider 1998). Instead, time delays 
are well measured for some quasars with great accuracy and 
arc less affected by spurious (and uncontrollable) effects. 
The unknown parameters to be determined are the source 
position {rs,9s), the shear quantities (7, S-y), the main lens 
model parameters (different for each model) , and the Hubble 
constant Hq. Actually, the higher complexity of the mod- 
els and the introduction of the time delays as constraints 
lead to more complicated equations than those considered 
in CCRPOl and CCRP02, and it is not possible to manipu- 
late them to reduce the set of lensing equations (2), (4) and 



(5) to a simpler form so as to speed up the calculations. Let 
us write equations (4) and (5) using the four images i, j, 
k and /, and the two equations coming from the time delay 
ratios among them, that are independent on h: 



A-f-obs ' 



Atu 



(16) 



where At°j, At"^ and At°i are the measured time delays. 
The introduction of these two other equations permits to rise 
the equations number, allowing to give a different constraint 
by the images positions that does not depend on the first 
derivative of the potential, but only on the potential. 

In summary, we consider the observables without the 
errors, using their mean values: while for the image positions 
this assumption seems immediately justified, for the time 
delays the uncertainties can be considerably large. In a next 
section, we will show more accurately which this assumption 
is reasonable. 

We have a number of 10 equations (8 from lens equa- 
tions and the 2 due to the time delay ratios), that we com- 
bine in a useful manner to reduce their number to the un- 
knowns number. For example, in the case of Model 1 the 
unknowns are 6, i.e., Ts, 9s, 7, 9^, b, and a, with in addi- 
tion the Hubble constant. We have to subtract the 8 lens 
equations to reduce their number to 4, and we then add the 
two time delay ratios to complete the system of 6 equations 
in 6 unknowns. We use the same procedure for each model 
to obtain a system of n equations and n unknowns. We can 
numerically solve this system of n equations to obtain n un- 
known parameters (for our models n is less than 10). Each 
result of the solution of the system will be a set of n values. 
When we go to solve the system, we do not obtain a single 
solution, but a set of solutions, since the system is highly 
non linear; these solutions have to be selected by means of 
some selection criteria in order to obtain values of parame- 
ters which give rise to physically plausible models. 

The search for the roots of the system begins with a 
choice of a range of parameter values: one must be sure that 
there are no roots outside the chosen range, and also a very 
large interval does not necessarily include all roots, increas- 
ing CPU time. We choose a physically acceptable interval 
for parameters, giving A/" random starting points to the al- 
gorithm, where A/" is a number fixed by the user^ . Prom these 
starting points the CPU tries to converge to the solutions, 
using the Newton's method. In particular, to assign the same 
probability to the starting values of a parameter in the rel- 
ative range, we generate these values uniformly distributed, 
also to avoid introducing any bias in the search. 

We have developed a simple code, named LePRet Lens- 
ing Parameters Reconstruction, written for Mathematica. 
Generally, solving the equation system yields M. < Af solu- 
tions, because for some values of the starting points New- 
ton's method fails to converge to a solution. 



^ Af should be large enough to explore a wide region in the pa- 
rameter space, but not too large so as to minimize CPU time. 
A possible strategy is to fix A" to a suitable value (for example, 
2000 or 10000) and run the code more than one time, 
t 'Lepre' is the Italian word for 'hare'. 
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4.2 Selection criteria 

These M solutions are not all physically acceptable and to 
sort among these we have to impose a set of selection crite- 
ria, that constrains parameter values to physically plausible 
ranges. Schematically, we can describe these selection crite- 
ria as following. 

(i) < T's < rnax{ri,rj,rk,'ri}: a lens will not produce 
images arbitrarily far away from the center of the lens; for 
large values of rs, there will be one image only at {r,9) = 
{rs,6s), and possibly others near the center of the lens; in 
particular, we impose this cut because if the source is outside 
the ring delineated by the most distant of the images it is 
not possible to generate a 4-images configuration^ . 

(ii) < 7 < 7cr«: the shear magnitude is positive by 
definition. For the separable models, we choose as upper 
limit 7 w 7crit, where "fcrit is defined such that for values 
of 7 > 7crit the estimated Ho becomes null, and depends 
on the particular lens system to be considered (Wucknitz 
2002). For PG 1115+080 is 7crit = 0.22, instead for RX 
J0911+0551 is "fcrit = 0.56. For axially symmetric lenses 
we cannot fix an upper limit (since the hypothesis made in 
(Wucknitz 2002) does not work), but in the analysis it is 
possible to obtain estimated value of 7 higher than the one 
obtained for an elliptical separable model. 

(iii) The shear is approximately well directed: we mean 
that the position angle d-y must be directed towards the ex- 
ternal mass disturbance; e.g., if the shear were due to an 
external group of galaxies, then 6^, should be aligned with 
the cluster mass centre, since there are no reasons why it 
should point elsewhere. For the axially symmetric models 
this cut is not necessary because an exact value of 9-y is 
selected automatically from the system solution to account 
for image configuration. Instead, for an elliptical model, we 
have to take into account the degeneracy existing between 
external asymmetry (the shear 7 and relative orientation 0^) 
and internal asymmetry provided by axial ratio q and angle 
Oq. Quantitatively, one could accept only values comprised 
in the range {6g min,dG max), where Oamin and Oomax are 
respectively two extreme galaxies that bound the external 
group of galaxies. 

(iv) The elliptical profile of the galaxy has to be physically 
plausible: we assume qo < q < 1, where qo is a particular 
value of q such that, for fixed value of a, the isodensity 
profile is elliptical or nearly elliptical, and not strongly boxy; 
so, we avoid those kinds of solutions to select the physically 
plausible one. 

(v) The range for other galaxy model parameters is chosen 
to have plausible values of k. For Model 1 and Model 3 we 
must have < a < 2; as a matter of fact, the surface mass 
density scales as r™~^, so that a < 2 is needed in order to 
have K monotonically decreasing with r. On the other hand, 
we consider a > in order for the projected mass inside r 
(that scales as r") to be reasonable. Since 1/) > 0, we have 
to assume the condition 6 > for the strength parameter. 



§ One can verify it using tlic web tool developed by K. Rat- 
natunga, which generates tlie images of a lens system once given 
the lensing potential, the source coordinates and the observational 
characteristics (see http://mds.phys.cmu.edu/ego_cgi.html). 



(vi) Plausible index ind(K). The index is the logarithmic 
derivative of k, i.e., ind{K) = . A lower bound can be 
fixed considering that k has to be higher than the luminous 
profile; an upper bound for PG 1115-1-080 is fixed following 
(Williams & Saha 2000), i.e., ind{K) < —0.5. For power law 
model this selection criterium reduces to a < 1.5. 

(vii) The set of parameters so found solves lens equations 
and time delay ratios equations. We insert the solutions in 
solved equations, to check if the solution found is the correct 
one or a result from wrong convergence of the numerical 
algorithm. We impose that the values of parameters so found 
verify the equations within an accuracy fixed by the user. 

(viii) hmin < h < hmax- From the three time delays we 

obtain three values of h (hij , hik and hu)*^ , and eliminate all 
the solutions such that the predicted values of these three 
h's differ more than an e {— 0.1) from each other. Then, we 
estimate h = {hij+hik+hu)/3, selecting only those solutions 
which give rise to values of h in the range {hmin, hmax), be- 
ing these two parameters fixed by the user. An upper bound 
is given by age of globular clusters and by the limits from 
nucleochronology, which indicate an age of the Universe of 
to > 7.8 Gyr; instead, a lower bound on h can be given 
by to 20 Gyr, since we do not observe stellar systems 
with ages greatly exceeding this value. In particular, to be 
conservative, we choose {hmin, hmax) = (0,2). 

Obviously, one can also change the order of the constraints, 

and remove or add some of these. At the end of the applica- 
tion of such selection criteria, and after having run the code 
more times, we get a set of C solutions. It is clear that there 
is a family of different combinations of the parameters, that 
is consistent with the observables and the selection criteria. 



4.3 Statistical interpretation 

It is possible to interpret the final sample of the solutions 
for each parameter using the concept of Bayesian proba- 
bility. In this particular framework, the lensing parameters 
and Ho represent the "variables" of the "system", that are 
linked by a series of "relationships" , that need not be single- 
valued, i.e., the lens mapping and time delay equations. We 
can associate to the "variables" a single value function, that 
is usually named "information", or more familiarly "prob- 
ability density function", which indicates how "likely" the 
particular combination of the parameters is. Our final sam- 
ple of solutions, collected into the histograms, determines all 
the amount of "information" that we obtain about the "sys- 
tem" itself. This sample is a subset of the parameter space, 
that contains all we have to know about the "system". To 
extract a final result, we have to select a better estimate of 
each parameter, making a reasonable choice. We could use 
a figure of merit as made by (Saha & Williams 1997), but 
here we use the mean value of all the collected values as final 
estimate for each parameter, instead for the uncertainty on 
this estimate we think to be conservative in the choice of the 
68% confidence level, i.e., the range in which the 68% of the 
solutions is included. In the first place, this choice works well 

^ We estimate three different h to give more freedom to the ob- 
tained solution and to take into account that we are using a nu- 
merical method. 
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in the simulations. Then, we argue that our choice allows to 
take into account the weight of each single result in the sam- 
ple. When the distribution of the values is symmetric, other 
choices as the maximum or the median are equivalent, but 
in some cases the histograms have a dogroo of asymmetry, 
different for each lensing model and quasar system, due to 
the particular configuration of the images, the degeneracies 
and the available data. 

As we said before, by means of a simulated system, we 
will sec that this estimate of the uncertainty in the parame- 
ter determinations in fact allows us to recover the values of 
these ones. 



5 A 'TWO-SYSTEM' MODEL 

In a recent paper, Saha & Williams extract an estimate of 
Ho by means of the joined use of more lenscd systems (Saha 

6 Williams 2004). Here, we show how it is possible to build 
a 'two-system' model fitting two lens systems with general 
lensing potentials and a shared Hq. In order to obtain a 
simpler solution to this problem we use an elliptical potential 
with a not fixed angular part Funk{6) and external shear 



1p = Funk{0)r" + Ipshear, 



(17) 



already used in CCRP02. As shown in that paper, if we do 
not assign the angular part of the potential, it is possible to 
write the time delay between two images i and j in a simple 
form and, in particular, for hij the relation 



h 



— 1 ''"100 



[(«-2)(r? 



■rj) 
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+{a - 2)7(r ■ cos 2{ei 
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holds, where the lens parameters related to the angular trend 
of the lens galaxy do not appear, but only the source posi- 
tion, the shear, the parameter a, and the normalized Hubble 
constant hij . The choice of this potential allows us not only 
to explore a wide range of models but also to obtain equa- 
tions with a more little number of unknowns and then more 
simply solvable. 

We impose that the h' s due to different pairs of 



images and the two different systems (i.e., h 



(1) 



h 



(1) 



h 



(1) 1,(2) 



13 ' ^ife'' ^iD" ^-re equal. We could solve a 



system of equations of the form h 



(1) 



^3 ' 

s , Cs ,7 



(2) 



in order to obtain the 10 parameters 



,eF,7^'',^f,a^"- We verified, 
anyway, that it is better to introduce a figure of merit of the 
form 

r = ih^ - h^r + ih^ - & + ih^ - h^f + (19) 

and we minimize it. We calculate the derivatives with re- 
spect to the 10 lens parameters and solve this system of 10 
equations. Following the previous section, we impose some 
selection criteria to select the physical solutions. In this case, 
in order to reject many unphysical solutions we have to im- 
pose a criterion not used previously: i.e., we require that the 



source has a more constrained position . This criterion is 
not arbitrary, since the source can be located only in a lit- 
tle region to be able to generate a particular configuration 
of the images (see (Saha & Williams 2003)). After having 
obtained a sample of solutions we have to calculate the six 
h's by different pairs of images and the two systems, and 
perform the mean of these values; finally, we impose the 
constraint on h's and obtain the final sample of solutions, 
that gives out the estimates of the parameters. 

We want to stress the fact that the algorithm discussed 
in this Sec. and hence the developed method, has as main 
issue the estimate of h. It does not allow to obtain precise 
estimate of the position of the source and of the orienta- 
tion of the external perturbation, since we impose a strong 
constraint on them, primary in order to obtain the correct 
estimate of h. As we will see later, it also gives us a good 
determination of the other parameters (i.e., 7, a, and h). 

We will compare the results with those obtained by us- 
ing a similar function of merit to model the single lens sys- 
tems with the same elliptical potentials. In this case the 
function F is simply modified introducing only the terms 
with the index (1). 



6 APPLICATION TO SIMULATED SYSTEMS 

In this section we verify the correct working of the codes, 
simulating real systems with imago positions and time de- 
lays exactly known (i.e., without observational uncertain- 
ties). Later on, analyzing the system PG 1115-1-080, we will 
show that this last assumption is reasonable, also for the 
real systems. 



6.1 Simulations 

For each simulation we fix the system parameters and obtain 
the image positions and the three time delays by solving the 
lens equations (4) and (5), and evaluating the time delay 
function. We use these quantities as input for the codes in 
order to verify if they are able to recover the correct values 
of the originally fixed parameters. 

In these simulation we adopt a flat homogeneous uni- 
verse with cosmological constant, fixing: 



{Qm, ^A, nk,h) = (0.3, 0.7, 0.0, 0.7) 
and: 

(zd,«.) = (0.31, 1.722), 



(20) 



(21) 



giving rioo = 33.37 days arcsec^^. 

In addition to the expected solution, that is recovered 
with high accuracy, the application of the codes generates 
other solutions, different by the first one, due to the exis- 
tence of the degeneracies among the lensing parameters, that 
we will describe in the next Section. The statistics on this 
solutions generates our estimate of each parameter. There- 
fore, the uncertainties of the parameters in the simulations 
and the real cases are not generated by a sort of propagation 
of the errors of the observables, but are instead the result of 
the degeneracies. 



The indices at the exponents specifies the lensed system. 



We impose the constraint on the angle 6s 
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Table 3. Simulated and estimated parameters, the angles are written in radiants. 
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The application of the method to simulated systems 
allows to verify how well it works in recovering the lens pa- 
rameters and the normalized Hubble constant h. We resume 
the results in Table 3, where we report the simulated values 
of the parameters and the estimated ones. 

Finally, we show the distributions of the values of h 
for the 5 models, adding the total distribution, obtained 
averaging the single distributions, in Figure 1. In particular, 
we can obtain a "marginalized" estimate of h averaging the 
distributions, taking into account all the recovered values in 
the distributions, and giving them finite weights: if Ni(h) 
indicates the distribution of the i-th model, the combined 
N{h) is given through a mean; this final distribution takes 
into account the difi^erent weights for the different models. 

We obtain the global result hest = 0.69±0.11. We may 
in fact consider this procedure a way to obtain a final es- 
timate of the Hubble constant unaffected by spurious un- 
certainties of the single codes, having taken into account 
different kinds of models and the whole parameters space. 
In this way, we can see that the final estimate of h is affected 
by a lower uncertainty. 



0.8 
0.6 



I \ 

I \ 



Legena 

■ Model 1 

Model 2 

■ Model 3 

■ Model 4 
Model 5 
All 




Figure 1. Distributions of the recovered values of tlie Hubble 
constant. N indicates a number of values normalized to the total 
number of results obtained for each model. This definition for 
N will be used for all the following continuous distributions. The 
distributions are obtained interpolating the recovered histograms. 
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Table 4. Simulated and estimated parameters for the model with 
a not fixed angular part [Two-systems model). 
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Table 5. Simulated and estimated parameters for the model with 
F{Q) unknown, fitted to a single system. 





simulated 


estimated 


' s 


0.25 


"•^^-0.18 




1.6 


1 70+"'2s 

-^•"^-0.47 




0.13 


^■1^-0.10 




3.00 


■^■^^-0.33 


„(1) 


1.1 


1 20+°'^^ 


h 


0.7 


"•""_0.18 



We have to verify whether the algorithm for the model 
with a common h is able or not to recover the lens param- 
eters and, in particular, the Hubble constant. The results 
are reported in Table 4: it is possible that some parameters 
are not perfectly recovered, but the main ones arc obtained 
with good accuracy. In addition, we verify the correct work- 
ing of the code when we fit the same potential to a single 
simulated system. The results are reported in Table 5. We 
note that the estimate of h for the two-system model has a 
lower uncertainty than that of the single modelling 

6.2 Model reconstruction from other lens models 

It is well known that the main uncertainty in applying the 
gravitational lensing in order to estimate the Hubble con- 
stant and to investigate the astrophysical properties of the 



galaxies comes from the lack of knowledge of the "true" 
galaxy model (see, for instance, (Schneider, Ehlers & Falco 
1992)). Here, we want briefly to analyze the biases and the 
errors due to the use of an "incorrect" parametric model to 
fit the observable data. A way to face this important ques- 
tion consists in to build simulated systems for a lens model 
and then to fit the image positions and time delays generated 
by it, using an other functional form, studying the change 
in the lens parameters and above all in h. We can furnish 
a qualitative estimate of the effect of the model dependence 
on Ho- At the same time the procedure allows to quantify 
this "systematic" errors. 

By means of a simulated system built using the Hub- 
ble model, wc fit to the image positions and the time delays 
so obtained the Model 1 and Model 2. The analysis of the 
most significant parameters shows interesting trends which 
are partially already well known. In particular, in our sim- 
ulated system h — 0.7, and the fitting of the other models 
gives us the means values h = 0.35 for the Model 1 and 
h = 0.46 for the Model 2, respectively with a percentage 
change of 50% and 34%. This percentage obviously changes 
if we simulate other lens systems. This shows that if the 
"correct" model for a lens is one with constant mass-to- 
light ratio and it tries to shape it with a separable model 
we obtain a lower estimate of h than that obtained with 
the first one; this verifies some previous results in literature 
(see (Kochanek 2002) for an analysis made on real systems). 
Similar trends are obtained if we create a simulated system 
using a de Vaucouleurs model. It is also possible to analyze 
the uncertainties introduced by the lack of the internal el- 
lipticity of the lens galaxy. If we try to fit with the Model 
1 the observables generated with the Model 2 we obtain a 
lower mean /; but if we consider the errors this estimate is in 
agreement with the simulated value, instead the estimated 
value for a raises. 



7 PARAMETER DEGENERACIES 

By means of simulated systems we can also obtain statis- 
tical correlations among parameters in order to investigate 
the degeneracies among them and to study the effect of vary- 
ing each parameter on the final Hubble constant estimate. 
Below, we discuss the results we found for each model. 

• Model 1. We see that rs,7,&, and h correlate each 
other, and anticorrelate with a. In particular, we verify the 
scaling laws (Wucknitz 2002): 

Ts 'X.2 — a, 70c2 — a , hoc 2 — a. (22) 

Wucknitz & Refsdal 2001 give a simple interpretation for 
these scaling laws in terms of mass-sheet degeneracy. In Fig- 
ure 2 we report as an example the correlation h — a, fitting 
the line h — a{2 — a), where a is a proportionality constant. 

• Model 2. For this isothermal model we observe that 
rs,b, and h correlate among each other, and anticorrelate 
with 7 and g; almost all of these are linear correlations with 
a high correlation degree. We note (see also later considera- 
tions about correlations for Model 3) that the absence of the 
parameter a changes the correlations among the shear 7 and 
the other parameters, i.e., since a is fixed, 7 now correlates 
negatively with Vs, b, and h, and not positively. 
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Figure 2. Scaling law /i oc 2 — a for Model 1. 

• Model 3. rs,^,q,b, and h correlate among each other 
and anticorrelate with a. Here, the positive correlation 
among 7 and the other parameters comes back: this allows 
us to confirm that, in Model 2, it is the condition a = 1 
that changes the correlations of the external shear and not 
the presence of an intrinsic cUipticity. We also note the ob- 
vious positive correlation between q and 7 (as for Model 2). 
Higher values of the internal ellipticity (i.e., lower values of 
q) require lower values of 7. 

• Model 4. For the de Vaucouleurs model we find similar 
linear correlations: r^, 7, b, and h correlate among each other 
and anticorrelate with i?e, which now assumes the role of a. 
For example, for higher values of Re, this model predicts 
lower values of Hubble constant. 

• Model 5. For the Hubble model we find the same cor- 
relations of the preceding model, simply replacing Re with 
the core radius s 

For real systems it is more difficult to obtain these cor- 
relations. Therefore, we do not discuss the results for them, 
and we prefer to resume the obtained dependencies in the 
following: 

• The parameters r^, &, and h always correlate positively 
among each other; i.e., more massive lenses require higher 
values of the radial position of the source and of the Hubble 
constant. The observations show that more luminous galax- 
ies give larger angular separations of the images, in agree- 
ment with our correlations: a larger h generates a larger Vs 
and, hence, a larger separation of the images. 

• The parameters that determine the radial profile of the 
lens galaxy {a, Re, or s) correlates negatively with the other 
ones. 

• Except for the Model 2, all the models reveal a positive 
correlation of the shear with rs, 6, and h. 



8 APPLICATION TO PG 1115-|-080 AND RX 
J0911-I-0551 

Having checked that the numerical codes indeed work cor- 
rectly recovering the values of the lensing parameters and 
Hubble constant, we apply them to real quadruply imaged 
systems. We need a four images system with a good astrom- 
etry of the lensing galaxy and image positions, and time 
delays measured with high accuracy. Then, we also need 



that there is a single galaxy acting as lens: for instance, 
witli our models we cannot study a gravitational lens like B 
1608-f-656, since there are two lensing galaxies; in this case, 
a more complex (two) lens model should be used. There 
arc only three systems that satisfy all these requirements: 
PG 1115-F080, RX J0911-F0551, and B 1422-^231. Here, 
we choose to apply the codes only to PG 1115-1-080 and 
RX J0911-I-0551, since B 1422-1-231 has time delays mea^ 
sured with a high uncertainty (see, for instance (Patnaik & 
Narasimha 2001)). Wc will also combine the results from 
the two systems to get a better estimate of the Hubble con- 
stant. We adopt a flat cosmology with (Clm, ^a) = (0.3, 0.7), 
discussing for each system the influence of changing the cos- 
mological parameters on the final estimate of Hq. 

8.1 Application to PG 1115-|-080 

PG 1115-1-080 was discovered by (Weymann et al. 1980), 
originally as a triple lensed quasar. Later, it has been possi- 
ble to split an image (the image A) in two components, Ai 
and A2. Therefore, this system consists of four images (^1, 
A2, B, and C) of a radio quiet QSO at Zs = 1.722, with 
an elliptical galaxy as lens belonging to a group of galax- 
ies at Zd = 0.31. This group, situated at South- West, con- 
tains ~ 11 galaxies, with a luminous centroid at {rg,6g) = 
(20"±0.2",-117''±3°). Iwamuro et al. fitted the luminous 
profile of the lensing galaxy with a de Vaucouleurs model 
with Re = 0".58±0".05, and measured an ellipticity ~ 0.1 
and a position angle of ~ 65° from north (Iwamuro et al. 
2000). 

Here, we use image coordinates measured by (Impey et 
al. 1998), and the time delay obtained by (Barkana 1997). 
In particular, in (Barkana 1997) it measured a time delay 
between the components B and C of Atsc = 25. Oil. 7 days, 
consistent with the previous value 23.7±3.4days from 
(Schechter et al. 1997). By contrast, the time delay ratio 
rABC = AtAc/AtBA = l-13lo!i7 fouud by Barkana is in 
conflict with the value 0.7±0.3 from (Schechter et al. 1997). 

We can fit this system with the models previously dis- 
cussed and observe some peculiar trends in the obtained 
parameter values, that we discuss in the following: 

• The source positions are consistent with each other ex- 
cept for Model 2, that has a little disagreement in its mean 
value. 

• As expected, axially symmetric models predict a higher 
mean value of the shear, to account for the lack of an internal 
asymmetry, being in agreement with the previous estimates 
of ~ 0.1. It is also interesting to note that 6-, is perfectly ori- 
ented towards the luminous centroid of the external group. 

• The isopotential profile has a small ellipticity, which 
increases when we go to consider the relative isodensity pro- 
files. In particular. Model 5 predicts an EO-El galaxy, while 
Model 3 is also consistent with more elliptical profiles. In- 
stead, 6q marginally agrees with the luminous profile orien- 
tation, but is nonetheless in agreement with results obtained 
with other techniques (Impey et al. 1998). 

• We obtain a values consistent (also if marginally) with 
the nearly isothermal model, with major preference for 
Model 3 (see (Kochanek 2002)). 

• Model 4 predicts a value of Re in agreement with the 
measured one of 0".58±0".05 by (Iwamuro et al. 2000), while 
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Table 6. Estimated parameters for PG 1115+080. 9q is oriented 
from West to North. 
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for Model 5 the core radius s is very small, being consistent 
with a null core -within the uncertainties. 

• At least, the most important result concerns the es- 
timated Hubble constants. Constant M/L models predict 
higher values than separable models, as yet found previously 
in literature (Courbin et al. 1997; Kochanek 2002). Then, for 
the isothermal Model 2, we verify a result of (Impey et al. 
1998) that gives a lo-w value using a SIE, an isothermal model 
similar to our model, but -with a different angular part. 

We report the results of the application of our pro- 
cedure in Fig. 3. In Fig. 4 we group together the differ- 
ent distributions. The mean distribution gives the value 
Ho = 58±27 kms"^ Mpc"^ The mean distribution for 
the constant mass-to-light ratio models gives us Hq = 
73±22 km s~^ Mpc~^ in according with the recent results 
in (Kochanek 2002) . If we change the values of cosmological 
parameters, our final estimate can change of ~ 3%, but the 
spread becomes ~ 12% if we also consider inhomogeneous 
models. 

Finally, in order to verify the honesty of our assumption 
of considering the observables without errors, we proceed 
as following: we consider as input values for the codes the 
mean values of the observables, and then, we change the 
time delays of an amount about the 15%. For example, for 
the Model 1 we obtain the results in figure 5. The two best 
values of h are: hdown ~ O.SOtg-ji and hup — 0.64to;25) 
in agreement with each other and with the value obtained 
using the central value. Therefore, the uncertainty in the 
time delays is extensively absorbed by the uncertainty due 
to the parameter degeneracies. 




Figure 4. Distributions of all the recovered values of the Hubble 
constant for PG 1115+080. 
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Figure 5. Histograms of the recovered values of the Hubble con- 
stant for PG 1115+080 with Model 1 changing the time delays. 



8.2 Application to RX J0911+0551 

RX J0911+0551 was discovered by (Bade et al. 1997) as 
a quadruple imaged QSO, in the ROSAT All-Sky Survey 
(RASS). The source QSO is situated at Zs = 2.8, while 
the lensing galaxy is at Zd = 0.77. The image configura- 
tion is peculiar: the three images A^, A2 and A-^ are close 
to each other and to the lensing galaxy (Ai — A2 = 0.478", 
A2 — A3 = 0.608"), while image B is more distant than the 
lens, requiring a large external shear. The only intrinsic el- 
lipticity does not allow to account for this configuration, and 
we need an external asymmetry. There is, in fact, a nearby 
group at about 38" South- West, and at redshift Zgroup = 
0.7±0.1 that can be the origin of this external perturbation. 
(Hjort et al. 2002) measured the time delay obtaining the 
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Figure 3. Histograms of the recovered values of the Hubble constant for PG 1115+080. The area under the histograms is normalized to 
unity. 



results AtsAi = — 143±6 days, At_BA2 = — 149±6 days and 
AtflAa = — 154±16 days. It is also possible to observe a sec- 
ond galaxy near the main one, that can afTect the modelling. 

In the following we present the main results of the ap- 
plication of our codes: 

• The models predict nearly the same values for the 
source positions, except Model 3, that provides a lower value. 

• The axially symmetric model yields a high value of the 
shear; for Model 2 we obtain 7 = 0.22, while Model 5 requires 
a lower value. The shear angle points towards the external 
group, allowing for the particular configuration observed. 

• The elliptical model requires a density distribution with 
high ellipticity consistent with an E1-E5 galaxy, with a po- 
sition angle oriented towards the external group. 

• Model 1 gives a lower value of a in agreement with the 
value from (Schechter 2000), while for the other elliptical 
model this value is high to account for the low values of 7 
and q. 

• The effective radius Re is low (~ 0.2), and the core 
radius of the Hubble model is consistent with zero. 

• The recovered value of Ho are higher than those ob- 
tained using PG 1115-1-080, with dramatically high uncer- 
tainties. The distributions of these values are too much flat 
giving us a little quantity of statistical information; maybe, 
using more complex models and other constraints we could 
improve these results. 

Our method allows to obtain the better estimate for 
the system parameters, but this circumstance does not also 
allow us to claim that a particular model fits the image 
positions and the other observables with higher accuracy. We 
note that the axially symmetric models are able to recover 
the correct position of the images A\, A2 and A^, but do not 
allow to get the position of the fourth image due to lacking 
of internal ellipticity. 

We may also try to apply an iterative procedure to ver- 
ify the correct working of the codes for this lens system. We 
fix the values of almost all the parameters leaving two of 
them free. Then, we choose one of those two, solving the 
lens equation relative to the image B and obtaining an es- 
timate for the other parameter^^. After finding this value, 
we can iterate the process solving the equation for the other 
parameter. After that, we again calculate the image position 
without being able to recover the correct ones. Therefore, we 



tt We solve this equation with 1 unknown, fixing as starting point 
the value obtained using the code. 



Table 7. Estimated parameters for RX J091 1-1-0551. 
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argue that is not possible to fit the image positions with an 
axial symmetric model with external shear. 

Our results agree with (Burud et al. 1998); they show 
that the particular image configuration requires a minimum 
ellipticity for the galaxy of e^in = 0.075 and a minimum 
shear of = 0.15, applying the method of (Witt & Mao 
1997). In particular, we predict very high values of 7, except 
for Model 2 and Model 3, that needs a lower mean value, but 
in agreement with that lower bound. 

The histograms of the Hq values are shown in Figs. 
6. We collect all the distributions in Fig. 7. We give 
two final estimates of Hubble constant, the first one only 
including the elliptical models, and the second one in- 
cluding all models. Using the elliptical models we obtain 
Ho = 81±41 kms~^Mpc~^, instead taking into account 
all the models we have Ho = 77±43 kms~^Mpc~^. Using 
a power-law model with external shear (the " Yardstick" 
model) (Schechter 2000) gets a low value of a, in agree- 
ment with the value obtained for Model 1, and Ho ~ 
42 kms"^Mpc"\ with an uncertainty ~ 10 - 20%. This 
value does not agree with the result obtained here be- 
cause in (Schechter 2000) is used a time delay of 200 days 
among the mean image A and the image B, different by 
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the one we use. Instead, including in the model the main 
lens galaxy, the cluster of galaxy, and individual galax- 
ies in the cluster (Hjort et al. 2002) get an estimate 
of Ho = 71±4 (random, 2(t)±8 (systematic) kms~^Mpc~^, 
that agrees with the value we obtain. 

Fixing other values for cosmological parameters, the 
maximum spread is ~ 10%, but considering a Dyer & Roeder 
universe we have a value ~ 30 — 40%, since the redshifts 
Zs and Zd of RX J0911+0551 are higher than those of PG 
1115+080. Actually, RX J0911+0551 is not the ideal system 
to obtain an accurate estimate of Ho, since there is a high 
uncertainty introduced by the arbitrary choice of different 
cosmological parameters. 

Finally, we remark that a simple elliptical potential with 
external shear, also being able to account for image configu- 
ration, is an approximate attempt to model a complex sys- 
tem as RX J0911+0551. Our modelling, in fact, allows to 
obtain useful information about the system, but we expect 
to obtain better results using, for example, a SIS or SIE 
to describe the contribution of the external group, and not 
an approximation as the external shear. Moreover, it is nec- 
essary to take into account the contribution of the second 
galaxy and the other objects in the group. 
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Figure 7. Distributions of all recovered values of the Hubble 
constant for RX J0911+0551. 



Table 8. Estimated parameters for the fitting of the elliptical 
potential with not fixed angular part (two-systems model), and 
comparison with the results obtained fitting the same potential 
to the single systems. 



8.3 Marginalized estimate of Ho 

For each systems we obtained an estimate of Ho performing 
a mean of the final distributions obtained fitting each model, 
and then deriving the mean value and the uncertainty from 
it. Now, let us combine these results^, multiplying the two 
final distributions. The final estimate turns out to be : 

Ho = 56±23 kms"^Mpc"\ (23) 

The uncertainties in each estimate and in the final one are 
high (also if reduced in this last one with respect to the 
single systems), since it has to take into account all the 
degeneracies; however, adding other physical constraints and 
strengthening some of those already used, it is possible to 
reduce the errors ulteriorly. Also for the lens parameters 
there is this kind of problem, above all in the orientation of 
the lens galaxy and the strength parameters of some models. 
We think that this kind of "marginalized" estimate of Hq 
over a large sample of models can help us to overcome the 
problem of the lack of a correct independent knowledge of 
the lens model. 
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8.4 Application of the two-system model 

In addition to the usual analysis done for each system, we 
apply to PG 1115+080 and RX J0911+0551 the new method 
to simultaneously fit the general elliptical model ip = F(6)r" 
with external shear. In Table 8 we collect the results, adding 
the values obtained for the single systems fitted with the 
same potential. The uncertainties in the estimated value 
for RX J0911+0551 are very high as well as for the other 
models fitted in the paper, showing the presence of diffi- 
culties in the fitting performed with these simple models. 
The angular trends for the two-system model and the single 

tt For RX J0911+0551 we only use the two elliptical models. 



systems arc similar and there is an agreement for the esti- 
mated as. The global estimate of the Hubble constant is 
Ho = 491?]^ Km s~^ Mpc~^. Its determination is mainly 
determined by the uncertainties in the estimate of Ho for 
PG 1115+080, since the distribution of the recovered for RX 
J0911+0551 is flat, giving us a little quantity of information 
(see, for instance. Fig. 8). 

These final values, but in particular the first one, 
are in good agreement with the previous ones in the lit- 
erature. The results obtained by some of us in previ- 
ous papers arc also consistent with the present ones: in 
CCRPOl it was obtained Ho = 56^11 kms"^Mpc"\ using 
an elliptical potential without external shear and only PG 
1115+080, while in CCRP02 it has been found a value of 
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Figure 6. Histograms of the recovered values of the Hubble constant for RX J0911+0551. 
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Figure 8. Distributions of all recovered values of the Hubble 
constant for the model withb F(9) unknown, applied to the two 
systems. 



Ho — 58lj^5 kms ^Mpc ^, fitting a general elliptical poten- 
tial to PG 1115+080 and B 1422+231. Using the pixel- 
lated lens method, Williams & Saha the results from PG 
1115+080 and B 1608+656 are combined to finally get 
H(j = 61±11 kms~^Mpc~^, in good agreement with our re- 
sult. Using a minimization applied to B1608+656, (Koop- 
mans & Fassnacht 1999) get Hq = 65lg kms~^Mpc~^, only 
marginally in agreement with our result. In (Treu & Koop- 
mans 2002) it is obtained Ho = 591^^+3 kms"^Mpc"\ 
modelling PG 1115+080 with two different components so 
as to describe the luminous part and the dark halo, and 
using information by stellar dynamics. 



9 CONCLUSIONS 

In this paper we have presented a numerical method able to 
estimate lensing parameters and Hubble constant for a wide 
class of models. The model parameters as well as the Hubble 
constant have been estimated using as constraints the image 
positions and the time delay ratios. We used two classes of 
models: separable elliptical models and constant mass - to - 
light ratio profiles, adding an external shear to take into ac- 
count the presence of an external group of galaxy. For these 
models we solved the system composed by the combinations 
of lens equations and two time delay ratios, selecting the so- 
lutions by means of suitable physical constraints. For each 
model and each parameter, we obtain an ensemble of values: 
we used the mean as better estimate and a confidence level 
of 68% as error. In order to reduce the uncertainty due to 
the lens models on the estimation of the Hubble constant, 
we marginalized over all the models collecting the complete 
data set and obtained a final estimate of Ho and its error. 



which does not result dramatically underestimated, in this 
way, because of an a priori choice of the model. To test the 
code we created simulated systems, being able to recover the 
correct values of parameters. 

After the encouraging testing of the codes, we have then 
applied them to two real systems for which a measure of time 
delays has been possible: PG 1115+080 and RX J0911+055. 
For PG 1115+080 it was possible to get an estimate of 
Ho = 58±27 km s~^Mpc~^, consistent with other results 
in the literature, obtained using different techniques. For 
example, (Courbin et al. 1997; Keeton & Kochanek 1997) 
show that the isothermal and pseudoisothermal models pre- 
dict low values of Ho , while the constant mass - to - light ratio 
ones generate higher values. We can verify these results us- 
ing Hubble and de Vaucouleurs models, finding that a simple 
elliptical isothermal model as Model 2 predicts a very low 
value of Ho, in agreement with the value ~ 40 km s~^Mpc~^ 
obtained in (Impey et al. 1998). For RX J0911+0551, only 
the elliptical profile allows to fit the image configuration. 
These trends are also confirmed by the simulations. 

As previously said, we "marginalize" over the models 
since we do not know the "correct" form of the lens model, 
and hence we have thought to overcome this difficulty in this 
way. The combination of the two final distribution can help 
to reduce the uncertainties and to obtain more information 
from more lens systems, thus should avoid the problems in 
the fitting of the single systems. The combined estimate is 
Ho = 56+23 km s^^Mpc"^ 

The uncertainty in the final estimate can be further re- 
duced adding other models and including in the statistics 
other lensed systems, consistently with the uncertainty ob- 
tained using other methods (see, for instance, (Williams & 
Saha 2000)). If we consider the contribution of the smooth- 
ness parameter a, the change in Ho can be very high and 
comparable with uncertainty in our estimates of Ho- For ex- 
ample, the variability for RX J0911+0551 is ~ 30 - 40% 
(using these particular cosmological models), with a compa- 
rable uncertainty in modelling. 

The general method developed in this paper can be used 
to do more, allowing to obtain an estimate of Ho that is as- 
sumed for hypothesis as the 'same' for all the lensed systems 
(see (Saha & Williams 2004)). We can in fact fit simultane- 
ously the two systems, using a general elliptical potential 
with a not fixed angular part and a shared Ho- The final es- 
timate is Ho ~ 49-11 Km s"'^ Mpc"'^, lower than the result 
obtained by means of the marginalization of the 5 models 
already analyzed in the paper, but in agreement within the 
uncertainties. The hypothesis of a common Ho is ambitious 
and very strong, since we don't know if different lens sys- 
tem can be fitted in the same manner by using the same 
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lens model and the same Ho- In fact, fitting different lens 
models, we have seen that different ones for the two lensed 
systems give us different values of Ho. But we think that a 
"more-system" model can help to obtain a reasonable esti- 
mate. 

Further improvements are possible. It will possible to 
use other different models, for which it is not possible to 
write the potential in a simple form, such as NFW pro- 
files ((Navarro, Prenk & White 1996; Navarro, Frenk & 

White 1997)), or more complex ones (also realizing more- 
system models). To take into account the different compo- 
nents of the lensing galaxy, we can use different profiles to 
describe their components: for example, it can be used a 
nearly isothermal model to describe the dark halo and a dc 
Vaucouleurs one to describe the luminous profile. Then, we 
could also use an exponential profile to account for a thin 
disk, that elliptical galaxies sometimes seem to have. It is 
necessary a more accurate modelling of RX ,10911+0551 in 
order to give a better bound on the estimated Ho for this 
system, since we checked that it furnish a little information 
and a little statistical weight. Then, it is possible to shape 
double lensed system, also using the flux ratio to better con- 
strain the lens models to systems: we have to use a function 
of merit similar to that used for the two-system model to 
allow to have a necessary number of equations. 

Finally, of course, the application of the single mod- 
elling (after marginalization) and other two-systems ones to 
other lens systems with measured time delays can allow to 
get a more and more accurate and precise estimate of Ho, 
eliminating biases and errors linked to the use of each lens 
model. 



ACKNOWLEDGEMENTS 

It is a pleasure to thank P. Scudellaro for the interesting 
discussions we had during the making of this work and the 
referee P. Saha for his significative suggestion in order to 
improve the paper. 

REFERENCES 

Bade, N., Siebert, J., Lopez, S., Voges, W., Reimers, D. 1997, 

A&A, 317, L13 
Barkana, R. 1997, ApJ, 489, 21 
Blandford, R., Kochanek, C.S. 1987 ApJ, 321, 658 
Blandford, R., Narayan, R.. 1986 ApJ, 310, 568 
Bradac, M., Schneider, P., Steinmetz, M., Lombardi, M., King, 

L.J., Porcas R. 2002, A&A, 388, 373 
Burud, I. et al. 1998, ApJ, 501, L5 

Caldwell, R.R., Dave, R., Steinhardt, P.J. 1998, Phys. Rev. Lett., 
80, 1582 

Cardone, V.F., Capozziello, S., Re, V., Piedipalumbo, E. 2001 

A&A, 379, 72 (CCRPOl) 
Cardone, V.F., Capozziello, S., Re, V., Piedipalumbo, E. 2002 

A&A, 382, 792 (CCRP02) 
Chang, K., Refsdal, S. 1979 Nature, 282, 561 

Courbin, F., Magain, P., Kceton, C.R., Kochanek, C.S., Vander- 
rist, C, Jaunsen, A.O., Hjorth, J. 1997 A&A, 324, LI 

de Vaucouleurs, G. 1948 Ann. D'Ap., 11, 247 

Demianski, M., de Ritis, R., Marino, A. A., Piedipalumbo, E. 2003 
A&A 411, 33 

Dyer, C.C., Roeder, R.C. 1972, ApJ, 174, L115 



Dyer, C.C., Roeder, R.C. 1973, ApJ, 180, L31 
Dyer, C.C., Roeder, R.C. 1974, ApJ, 189, 167 
Hjorth, J. et al. 2002, ApJ, 572, Lll. 

Impey, CD., Falco, E.E., Kochaneck, C.S., Lehar, J., Mcleod, 
B.A., Rix, H.-W., Peng, C.Y., Keeton, C.R. 1998 ApJ, 509, 
551. 

Iwamuro, F. et al. 2000, Publ. Astron. Soc. Japan, 52, 25. 

Kassiola, A., Kovncr, L 1995 MNRAS, 272, 363. 

Keeton, C.R., Kochanek, C.S. 1997, ApJ, 487, 42. 

Keeton, C.R. 2001, astro-ph/0102341. 

Kochanek, C.S. 2002 astro-ph/0204043. 

Kochanek, C.S., Schechter, P.L. 2003 astro-ph/0306040. 

Koopmans, L.V.E., Fassnacht, CD. 1999 ApJ, 527, 513. 

Koopmans, L.V.E., de Bruyn, A.G. 2000 A&A, 358, 793. 

Kovner, L 1987, Nature, 325, 507. 

Mao, S., Schneider, P. 1998, MNRAS, 295, 587. 

Narayan, R., Bartelmann, M. Lectures on gravitational lensing, 

in Formation of structures in the Universe, ed. Dekel, A., Os- 

triker, J. P.. 

Navarro, J.F., Frenk, C.S., White, S.D.M. 1996 ApJ, 462, 563. 
Navarro, J.F., Frenk, C.S., White, S.D.M. 1997, ApJ, 490, 493. 
Patnaik, A.R., Narasimha, D. 2001, MNRAS, 326, 1403. 
Peebles, P.J., Ratra, B. 2002 Rev. Mod. Phys., 75, 559P. 
Refsdal, S. 1964, MNRAS, 128, 295. 
Refsdal, S. 1964, MNRAS, 128, 307. 
Refsdal, S. 1966, MNRAS, 132, 101. 

Romanowsky A. J. et al.. Science Express Reports, 28 August 

2003 

Saha, P., WiUiams, L.L.R. 1997, MNRAS, 292, 148. 

Saha, P., WiUiams, L.L.R. 2003, AJ, 125, 2769. 

Saha, P., WiUiams, L.L.R. 2004, astro-ph/0402135. 

Schechter, P.L., Bailyn, CD., Barr, R., et al. 1997, ApJ, 475, L85. 

Schechter, P.L. 2000, astro-ph/0009048. 

Schneider, P., Ehlers, J., Falco, E.E. 1992, Gravitational lenses, 
Springer- Verlag, Berlin. 

Treu, T., Koopmans, L.V.E. 2002, MNRAS, 37, L6. 

Walsh, D., Carswell, R.F., Weymann, R.J. 1979 Nature, 279, 381. 

Williams, L.L.R., Saha, P. 2000, AJ, 119, 439. 

Weymann, R.J., Latham, D., Angel, J.R.P., Green, R.F., Liebert, 
J.W., Turnshek, D.A., Turnshek, D.E., Tyson, J.A. 1980, Na- 
ture, 285, 641. 

Witt, H.J., Mao, S. 1997, MNRAS, 291, 211. 

Witt, H.J., Mao, S., Keeton, C.R. 2000, ApJ, 544, 98. 

Wucknitz O., Refsdal, S. in Brainerd, T.G., Kochanek, C.S., 
eds., ASP Conf. Ser. Vol. 237 Gravitational Lensing: Recent 
Progress and Future Goals. Astron. Soc. Pac, San Francisco, 
p. 157. 

Wucknitz, O. 2002, MNRAS, 332, 951. 

Zhao, H.S., Pronk, D. 2001, MNRAS, 230, 401. 

http : //cf a-www . harvard . edu/cast les/ 



